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Abstract 

We recently have presented first physical predictions of a spatially hybrid model that follows the 
evolution of a negative streamer discharge in full three spatial dimensions; our spatially hybrid 
model couples a particle model in the high field region ahead of the streamer with a fluid model 
in the streamer interior where electron densities are high and fields are low. Therefore the model 
is computationally efficient, while it also follows the dynamics of single electrons including 
their possible run-away. Here we describe the technical details of our computations, and present 
the next step in a systematic development of the simulation code. First, new sets of transport 
coefficients and reaction rates are obtained from particle swarm simulations in air, nitrogen, 
oxygen and argon. These coefficients are implemented in an extended fluid model to make the 
fluid approximation as consistent as possible with the particle model, and to avoid discontinuities 
at the interface between fluid and particle regions. Then two splitting methods are introduced and 
compared for the location and motion of the fluid-particle-interface in three spatial dimensions. 
Finally, we present first results of the 3D spatially hybrid model for a negative streamer in air. 
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1. Introduction 

1.1. Streamer discharges and single electron dynamics 

Streamers are ionized fingers that penetrate into a previous ly non- ionized region under the ac- 
tion of an electric field; they pave the way of sparks and lightning, but they also appear without 
these subsequent stages in corona reactors or in the form of huge sprites high above thunder- 
clouds; for recent reviews on streamer physics we refer, e.g., to [1, 2]. Streamers are character- 
ized by a thin charge layer around their head, which screens the electric field inside the streamer 
body and enhances it at the tip. This local enhancement of the electric field allows them to 
generate additional plasma ahead of the already ionized channel through impact of accelerated 
electrons onto neutral atoms and molecules. 
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As streamers in air at standard temperature and pressure contain at least of the order of 10 7 
electrons already at the moment when they emerge from avalanches through the build-up of a 
space charge layer [3, 4], it is extremely challenging to impossible to follow the growing number 
of electrons individually inside the streamer discharge during the further evolution. Therefore 
most streamer simulations are performed in a density approximation for the electron and ion 
dynamics in so-called fluid models. As reviewed below and in another paper in this issue [5], 
the numerical solution of fluid streamer models is very challenging as well due to the widely 
differing length scales of the problem; these simulations give realistic results for the macroscopic 
dynamics. 

But there are certain physical phenomena that are difficult to impossible to simulate with a 
fluid model. These are related to either small electron numbers or to high electron energies. 

1.1.1. Small electron densities and density fluctuations 

Regions with very low electron densities can dominate the macroscopic discharge evolution 
in two cases, namely during the initial avalanche regime, and when streamer propagation is un- 
stable and susceptible to branching. In these cases, the presence or absence of single electrons 
can be visible in the macroscopic evolution. A single electron can start a discharge; and the 
stochastic distribution of discharge inception times at low voltages can probably be traced back 
to the stochastic distribution of single electrons during the inception phase. Furthermore, recent 
streamer experiments in very pure gases gave the first experimental evidence of avalanches cre- 
ated by single electrons [6, 7]; these avalanches do not necessarily lead to branching, but rather 
can give the discharge a feather-like structure with many parallel "hairs" on the main channel. 
(We remark in passing that avalanches are neither necessary nor sufficient for streamer branch- 
ing, as discussed in more detail in [1, 8], but they can accelerate the branching of an already 
unstable streamer head [5]. 

7.7.2. High electron energies and electron run-away 

The electron energy distribution in the ionization front can develop a long tail at high ener- 
gies that is not appropriately modeled by a fluid model, even if the fluid model is extended by 
additional terms as in [9]. In particular, individual electrons can run away [10, 11, 12]. And when 
the electric field exceeds a threshold, the majority of electrons will run away and the density or 
fluid model will break down completely [13, 14, 15]. 

Much recent interest in extreme electron energies comes from the observation of Terrestrial 
Gamma-Ray Flashes. These flashes were first detected by the Compton Gamma Ray Observatory 
satellite in 1994, and soon a correlation with lightning discharges was found [16]. This energetic 
radiation is now generally believed to be the bremsstrahlung of electrons with extremely high 
energies. Meanwhile, even a considerable content of positrons is found in these flashes [17]. 
But how the electrons gain these energies in the discharge, is under debate, candidates are either 
relativistic run-away electron avalanches or run-away electrons emitted from streamers, leaders 
or other lightning discharge processes. Hard radiation was also observed from rocket triggered 
lightning near the ground [18, 19], and from long sparks generated in laboratory with Mega Volt 
pulses [20, 21, 22] as well from streamer coronas generated by a pulsed voltage of 85 kV [23]. 
In the laboratory experiments and in the approaching lightning leader, it is clear that local field 
enhancement at electrodes and/or streamer tips is responsible for high electron energies. The 
same mechanism could be at work in Terrestrial Gamma-Ray Flashes. 
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1.1.3. Why to develop a spatially hybrid model 

Discharge inception, streamer feathers and streamer branching on the one hand and the run- 
away and continued acceleration of electrons in streamer discharges on the other hand give plenty 
of motivation to study the dynamics of single electrons in a streamer. But — as remarked above 
and further elaborated in the next subsection — neither particle nor fluid model are presently 
able to deal with these challenges. 

In this paper therefore the next step is taken in the development of a spatially hybrid simula- 
tion tool. The hybrid simulation treats the majority of electrons inside the streamer in a density 
approximation, which is appropriate and efficient as the densities are high and the local field is 
low; the fewer electrons in regions with lower densities and high field are treated in a full particle 
model. New steps in the present paper concern the consistent density approximation of the par- 
ticle model, and the construction and motion of the interface between particle and fluid region in 
full three dimensions. This model will allow to calculate how the dynamics of single electrons 
in the high field region influences the streamer motion, and how this high field region in turn is 
generated by the growing streamer body. 

1.2. The state of particle and fluid models of streamers 
1.2.1. Particle models 

Most particle simulations of plasmas in general, and of streamers in particular, are carried out 
with a particle-in-cell/Monte Carlo collision (PIC/MCC) model. The PIC/MCC model follows 
the free flight of single electrons and treats all the (important) collisions stochastically [24, 25, 
10, 26]. 

An important constraint of particle models is the number of electrons they can treat; the 
electron number can easily exceed limitations of computer memory and computational power. 
Particle simulations with real particles can only be performed for the early phase of streamers, 
such as the avalanche phase or the avalanche-to-streamer transition. To follow the streamers in 
later phases, super-particles are typically used, where one computational particle stands for the 
mass and charge of many real particles. However super-particles have their own drawbacks, not 
only through the lower resolution, but more importantly, because they cause numerical heating 
and stochastic errors, as shown in [27, 4], which can have a detrimental influence on the simula- 
tion results. 

Due to the heavy computational cost for a large amount of particles as well as for solving 
the 3D Poisson equation, a full 3D particle model of a streamer is difficult or in many situations 
impossible to realize. A compromising approach is to keep the positions and velocities of the 
particles defined in 3D space as in the standard PIC/MCC procedure, while calculating the elec- 
tric field in ID [28, 24] or 2D [26, 25], and then to describe the electron motion in the electric 
field by interpolation. Runaway electrons generated by streamers have been addressed in the 
work of Moss et al. [10] where the 3D electric field is simplified into a two-step function in ID: 
the high field value stands for the field at the streamer head and a low field value stands for the 
field far ahead of the streamer head. Chanrion et al. [12] recently published runaway simulation 
results in a 2D geometry with radial symmetry. Superparticles are used in both simulations, but 
an energy-dependent re-sampling technique for the superparticles is developed in [12], and a 
low weight can be attributed to energetic electrons which allows to study electron runaway in 
negative streamers with higher precision. 

In our particle model, that constitutes a part of our hybrid model, both the electric field and 
the electrons are defined and calculated in 3D. 



3 



7.2.2. Fluid models 

The fluid model approximates electrons by continuous densities and is therefore computa- 
tionally much more efficient. The fluid model is in most cases derived by averaging over the 
Boltzmann equation [29, 30, 31]. The transport coefficients such as mobilities and diffusion 
rates, and the reaction rates are approximated as functions of the local electric field or the local 
mean electron energy. (More details on the local field approximation and on the energy equation 
can be found in [32, 9].) 

3D fluid models have become available only recently [33, 34]; progress was from ID [35] 
and 1.5D [36] density simulations to various 2D [37, 38, 39] fluid descriptions developed in 
the last 20 years. While in 2D, the radial nonuniformity of particle densities and the local field 
enhancement are included in the fluid description, the streamer evolution after branching or the 
interaction of streamers with another [40] or with lateral walls requires a model in full 3D. 

Most streamer simulations are carried out with fluid models [37, 38, 41, 42, 39, 43, 44, 45, 46, 
47, 3, 48, 49, 7]. Solving the fluid equations numerically is not an easy task due to the multiscale 
nature of streamers. New simulation techniques have been developed [3, 48, 50, 51], a further 
review can be found [5]. A high resolution for the inner structure of the streamer head is needed, 
where the electron and ion densities decay sharply and where the electric field changes rapidly 
in space and time. These numerical challenges have been met by adaptive grid refinement [3, 50, 
40] . Improvement of the numerical techniques is only one aspect of the recent development in the 
fluid models; more physics is included by using more complicated and realistic plasma-chemical 
models [52, 53, 54, 48], better techniques of modeling electrode geometries [39, 48, 49] and 
efficient descriptions of non-local photo-ionization sources [55, 56, 47, 57, 58, 59]. 

1.3. Developing the hybrid model 

To combine the computational efficiency of the fluid model with the detailed physical de- 
scription of the particle model, we here take the next steps in the development of a model that 
is hybrid in space, coupling a particle description of the electrons in the region of high field and 
low electron density with a fluid description in the region with low field and high densities. We 
refer to [60, 11,9] for the detailed concept of the spatially hybrid model. 

In previous work, we have identified deviations between particle and naive fluid descriptions 
of planar streamer ionization fronts [60], and we have described how to construct the numerical 
interface between a particle and a fluid model in ID, and how to perform the hybrid simulation 
of a ID streamer front [61, 9]. In [1 1] we have presented first results of hybrid calculations in full 
3D, emphasizing on the acceleration of electrons to energies above 3.5 keV in a streamer head. 

In this paper, we focus on two aspects in the model development, namely on the calculation 
of transport and reaction coefficients for the fluid model from the particle model, and on the 
construction of the interface between particle and fluid model in a 3D simulation. 

First, a set of transport coefficients and reaction rates is derived from the particle model for 
the fluid model; this is done for air in section 2.1.3, and for pure nitrogen, oxygen and argon 
in the appendix. In our previous papers [60, 9], we already had calculated these coefficients for 
pure nitrogen, but there we fixed the elastic total cross section rather than the elastic momentum 
transfer cross section. Therefore the calculated coefficients were inconsistent with each other 
when different distribution for the scattering angles were used in the electron-neutral collisions. 

Second, two splitting methods are developed in section 3.2 and 3.3 to determine the inter- 
face between the particle and the fluid region. The first one, the column based splitting, is a 3D 
extension of the splitting method developed for planar fronts [9]. The first 3D results of this 
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method were reported in [11], here we give the numerical details. The second splitting method, 
the full 3D splitting, is an improvement of the first method. It generates model interfaces wher- 
ever needed, e.g., at both ends of a double headed streamer or even when avalanches or streamers 
break up in several parts. 

We also present new results for the 3D fluid model as well as for the 3D hybrid model in air. 

1.4. Organization of the paper 

Particle and fluid model are the basic components of the 3D hybrid model. We first give the 
details of numerical methods and their implementation in Section 2.1 for the particle model and 
in Section 2.2 for the fluid model. The coupling of the two models is discussed in Section 3 in 
which the general coupling procedure is given in Section 3.1, the methods to split the simulation 
domain into particle regions and fluid regions are then discussed in Section 3.2 and Section 3.3, 
and the definition of the buffer region is discussed in Section 3.4. In Section 4, we present the 3D 
hybrid simulation result for a negative streamer in air. We finish with some concluding remarks 
and suggestions for further research in Section 5. The transport coefficients and reaction rates 
for nitrogen, oxygen, and argon are given in Appendix Appendix A. 

2. Particle and fluid model 

We here lay the basis for constructing the hybrid model in 3D in section 3. In subsection A, 
we recall the essentials of the particle model with Monte Carlo procedure and input data for the 
differential cross sections; and we explain how we calculate transport and reaction coefficients 
for the fluid model in the most consistent manner from the particle model. We recall that the 
fluid model needs an extension in high fields to fit the particle model [9]. In subsection B, we 
discuss the fluid model and its numerical implementation and how to solve the Poisson equation. 
Finally, we present first simulation results of the fluid model in 3D on a Cartesian grid. 

2.1. Particle model 

2.1.1. The Monte Carlo procedure 

The particle model follows the standard procedure of Monte Carlo (MC) particle simulations 
in plasma physics; it is summarized here. The particle model tracks the motion of all electrons 
while treating the neutral atoms and molecules as a random background. Since the mobility 
of ions is two orders of magnitude smaller than that of electrons, ions are treated as immobile. 
As the ionization density in streamers is low at standard temperature and at standard pressure 
or below [2], only collisions with neutrals need to be taken into account. During a collision 
with a neutral, the electrons change velocity and energy (in the case of an elastic collision only 
velocity); the collisions are modeled as instantaneous. Between collisions, electrons follow the 
equation of motion x = qE/m; numerically this motion is solved with the leapfrog method 

x„ +1 = x„+Afv„ + i, (1) 

V„ + . = V„_i + Af-E(x„,f„), (2) 

2 2 m 

as illustrated in Fig. 1. Here q and m are elementary charge and electron mass, x„ = {x,y,z) n is 
the electron position at time f„, and v n+ i = (v. v , v y , v z ), 1+ i is the electron velocity at time t n+ i. 

The collisions of electrons with the neutral background are characterized by probability dis- 
tributions for the collision time, for the type of collision (elastic, inelastic and ionizing processes) 
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Figure 1: The schematic illustration of the Leapfrog scheme. 




Figure 2: Collision frequencies v(e) in air at standard temperature and pressure as a function of electron energy e. 
The frequencies of different collision types are plotted on top of each other such that they add up to the total collision 
frequency. From below, the collision frequencies are successively elastic, rotational and vibrational, electron exciting 
and ionizing collisions with nitrogen, elastic, attachment, rotational and vibrational, exciting and ionizing collisions with 
oxygen, and elastic, exciting and ionizing collisions with argon. When the collision frequency is smaller than maximal 
v(e) < y„ MV , artificial "null collisions" without physical effect are added to make the Monte Carlo procedure more 
efficient. 



and for the scattering angle, and in the case of an ionizing collision also for the distribution of 
energy between the two outcoming electrons. The actual collision events are sampled from the 
probability distributions through random numbers in a Monte Carlo procedure. The probability 
distributions are determined by the differential cross sections for the respective gas composition 
and by the gas density. The cross sections used in the present paper are summarized in subsec- 
tion 2.1.2. 

The cross sections and therefore the collision frequencies in general depend on the electron 
energy. But if the next collision time is calculated for each electron individually depending on 
its present energy, and if the change of energy during the flight is not taken into account, the 
calculation is both numerically expensive and inaccurate. Therefore the "null collision" [62] is 
used in most MC codes, which is actually a pseudo collision because nothing changes in this col- 
lision process. With the null collision technique, the next collision time is sampled by drawing a 
random number right after the previous collision; the probability distribution is based on the max- 
imum of the collision frequency v max over all reasonably occurring electron energies. Then when 
the collision time is reached, the actual energy e c of the incident electron is calculated. Another 
random number is drawn to determine the type of collision process (elastic, inelastic, ionizing 
or null) for an electron with e c . The null collision accounts for the probability that no collision 
occurs at the previously determined collision time, because the actual collision frequency for 
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electrons with energy e c is smaller than v max . The loss of electron energy is determined by the 
type of collision process, and the scattering angle is determined by the incident electron energy 
and sampled with another random number. If the collision is ionizing, another random number is 
drawn to determine the energy distribution between the two out-coming electrons. The scatter- 
ing angles of the out-coming electrons are determined by their energies as described in the next 
subsection. The electrons then follow their Newtonian (or eventually relativistic) trajectory up 
the next collision. 

Fig. 2 shows the total collision frequencies as a function of electron energy in air (modeled as 
78.12% nitrogen, 20.946% oxygen, and 0.934% argon) at standard temperature and pressure; the 
frequencies of elastic, rotational and vibrational, different exciting and finally ionizing collisions 
are plotted on top of each other such that they add up to a total collision frequency. This total 
collision frequency reaches a first peak around e = 2 eV due to the vibrational collisions. For 
growing e > 2 eV, it then first drops and then increases again up to a maximum v max at e = 150 
to 200 eV, after which it decreases. When the collision frequency decreases, the mean free path 
length of the electrons increases. If the electric field is so high, that the mean energy gain between 
collisions exceeds the mean energy loss during collisions, the electrons can continuously gain 
energy and run away. Electron run-away is one of the effects we want to track with our model. 

An important difference with most other PIC/MCC codes is that we use real particles instead 
super-particles. Since our particle model is designed for the region with high field enhancement 
where the electron density is relatively low, we intend to use our particle model with real parti- 
cles. Therefore, as long as one particle is identified with one real electron, the common problems 
of a PIC/MCC procedure such as numerical heating, stochastic errors and low resolution at es- 
sential regions cannot occur. 

The electric field is calculated using a fast Poisson solver, see more discussion in Sec- 
tion. 2.2.3. 

2.1.2. Data for differential cross sections 

The cross section data can be obtained through either one or more of the following means: 
(1) measurement from single-scattering beam experiments [63], (2) ab initio quantum theoreti- 
cal calculations [64, 65] or (3) inversion of swarm experimental data [66, 67]. Different cross 
section databases have been used by different authors, depending also on the stage of streamer 
development, for example, Dowds et al. uses the National Institute of Standards and Technol- 
ogy (NIST) electron impact cross section [68] in a particle avalanche model [25], Babich et al. 
uses the Evaluated Electron Data Library (EEDL) [69] for the simulation of high energy electron 
avalanches [70], and Li et al. and Chanrion et al. use the Siglo database in the streamer simula- 
tions [4, 26]. None of the cross section databases mentioned above covers the energy range from 
low energy (less than 1 eV) to very high energy (above 1 MeV). 

In our model, the cross sections for electrons with energy up to 1 keV colliding with nitrogen, 
oxygen or argon particles are taken from the siglo database [71], which includes all important 
collision processes in the streamer development. Above 1 keV, the Born approximation [72] 
is used for elastic collisions, a fit formula in [73] is implemented for the electronically excit- 
ing collisions and the Born-Bethe approximation [74, 75] is used for ionizing collisions; these 
approximations continuously extend the siglo database. The distribution of scattering angles is 
discussed further below. 

In an ionizing collision, the incident electron loses energy to liberate an electron from the 
neutral molecule or atom. There are several measurements of the angular and energy distribution, 
the so called doubly differential cross sections for ionization. The energy distribution of the 
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Figure 3: Normalized probability that the primary out-coming electron carries a fraction r= e p j(e - e,„„) of the total 
available energy. The incoming electrons have energies e = 10 (solid), 50 (dotted), 100 (dashed), 500 (dash-dotted) eV; 
ffon = 15.6 eV is the ionization energy of nitrogen. The probabilities are normalized over 100 Ar from r=0 to 1. 



secondary electrons ejected upon ionizing collisions also has a clear influence on the simulated 
particle swarms or streamer avalanches [76], since clearly electron run-away is more likely if 
one electrons carries away most of the energy after the collision while the other one is slow. The 
measurements of the secondary electron spectra by Opal et al. [77, 78] have become a standard, 
and their empirical fitting equation has been widely used in different simulation groups. Opal's 
empirical fitting is used in our calculation as well. 

In Fig. 3, we show the probability distribution cr(e p , e)/<r(e) that the primary out-coming 
electron carries the fraction r= e p /(e - e !0 „) of the available energy, for given incident electron 
energies e = 10 (solid), 50 (dotted), 100 (dashed), 500 (dash-dotted) eV; here e p is the energy 
of the primary electron (that takes more energy than the secondary electron), and e, „ is the 
ionization threshold energy of nitrogen. cr(e) is the total ionization cross section for electrons 
with incident energy e, and cr(e p , e) is the cross section for incident energy e and out-coming 
primary electron energy e p . The probabilities are normalized over 100 Ar for r e (0, 1). The plot 
shows that an incident electron with high energy is likely to keep most of its energy. 

The scattering angles of primary and secondary electron are determined by their energies [62] 

with 

cos^, = y/ei/(e - e ion ), (3) 

where i = 1,2 indicates the primary or secondary electron. If the primary electron takes most of 
the energy, i.e., cos^-,- * 1, it leaves the collision mostly in a forward direction; the secondary 
electron with low energy co&Xi « then flies in a direction perpendicular to the primary one. 

For the distribution of scattering angles in elastic collisions, analytic approximations are 
used as well in our model due to the relative lack of experimental data. Okhrimovskyy et al. 
has derived a formula from the integrated and momentum transfer elastic cross section data from 
Phelps and Pitchford [66] for nitrogen; it has been implemented in our model for the electron-N2 
elastic collisions ^ 

/(e '* } = 1 ~f n > (4) 

An [1 -^(e)cos^] 2 

where I(e,x) is the normalized differential cross section, e stands for the electron energy in eV 
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and x is the scattering angle with respect to the direction of the incident electron. An empirical 
fit of the parameter £ for nitrogen is 

0.065 e + 0.26 VI 12 VI 

£0) = p p, (5) 

1 + 0.05 e + 0.2 VI 1 + 40 VI 

where e = e/eV is the electron energy in dimensionless units. The formula derived by Surendra 
et al. [79] has been used in a particle simulation for electron-C>2 and electron- Ar collisions in [79, 
80], and is also used here in our model for electron-Ch and electron-argon elastic collisions, 

I(e,x) = T- • ( 6 ) 

47T [1 + e sin 2 (y/2)] ln(l + e) 

The experiment has shown that the electron scattering angles are different for different colli- 
sion processes [63]. A complete particle model could apply different differential cross sections 
for different excitational collision process. However, the differential cross sections are not avail- 
able for all essential collision processes in the literature, and the available data are mostly mea- 
sured for electrons at a few discrete energies and not in the form ready to be implemented in the 
particle model. Therefore, the electron scattering from the rotational, vibrational and electronic 
excitations is assumed to have the same scattering probabilities as from the elastic collisions. 
In an excitational collision, the electron scattering angle is first sampled based on the incident 
electron energy, then the energy loss in the excitation collision is subtracted from this energy. 

A methodological difference with our previous swarm calculation for nitrogen [60] should 
be noted. In [60] we used the elastic momentum transfer cross section 

cr m (e)=2n\ (1 - cos^f) l{e,x) s^X d X< ( 7 ) 
Jo 

from the siglo database, but we fixed it as the total elastic cross section. This lead to inconsis- 
tent transport and reaction coefficients between models with isotropic and anisotropic scattering. 
Now the total elastic cross section cr, 

cr,{e) = In I I(e,x) sin;^ (8) 
Jo 

is adjusted depending on the applied scattering model, such that a constant momentum and en- 
ergy transfer rate is obtained independently of the angular scattering model. The change has 
negligible influence on the electron transport and reactions for low energy electrons, but influ- 
ences become significant for electrons with higher energies, as shown in [81, 79] and also found 
in our study. 

When anisotropic scattering is applied in the particle model, the ratio of elastic momentum 
transfer cross section cr m and total elastic cross section cr, can be obtained from Eq. (7,8). The 
total elastic cross section can then be calculated from the elastic momentum transfer cross section 
with these ratios. The ratio for nitrogen is written as [82] 

cr,(e) 2£(e) 2 \ 1 - £(e) / 

where £ is the same as Eq. 5, and for oxygen and argon the ratio is written as [79, 80] 

2is W = _2_( 1 _MiiO\ (10) 

<t,{6) \n(\ + e)\ e )' 

where the function I(e,x) are from Eq. (4) and (6) respectively. 
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2.1.3. Calculation of transport and reaction coefficients in air 

The electron transport coefficients and reaction rates are calculated in particle swarm exper- 
iments and used in the fluid model to keep the two models consistent. Two sets of transport 
coefficients and reaction rates have been calculated by particle swarm experiments: flux coeffi- 
cients and bulk coefficients [83]. (The use of flux or bulk coefficients is discussed in the next 
section.) In a constant field E, the bulk coefficients are calculated from the swarm using the 
following equations: 

H(E)\E\ = , (11) 

'2 — <1 

1 In N e (t 2 )- ^ N e (h) 

a(E) = , (12) 

fi(E)\E\ t 2 -h 



D(£) = <rv>-<r><v>, 



(13) 



and for the flux coefficients: 



H*(E)\E\ = <v : > (14) 

, (E) _J ]aN e (t 2 )-\aN e (t 1 ) 

fi*(E)\E\ t 2 -h 
nv ™ «r(f 2 ) 2 ) - <r(f 2 )) 2 ) - «r( fl ) 2 ) - <r( fl )) 2 ) 

D(£) = m^) ' (16) 

= m-w = m-w t (17) 

a*(E)u*(E) a(E)n(E) 

where E is applied along z axis, N e (t) is the total number of electrons at time t, r = (x, y, z) and 
v = (y x , v y , v z ) are the position and velocity of electrons and (■ ■ ■ ) denotes the average over all 
particles. k\ is the nonlocal parameter in the extended fluid model (see Section. 2.2 or [9] for 
the definition of k\). The particle swarm experiments have been done for a range of electric field 
from E — 5 to 250 kV/cm for air at standard temperature and pressure. By using similarity laws, 
they can be easily converted to other temperature and pressures. 

The flux coefficients represent the real motions of electrons, for example, the flux mobility 
times electric field fiE represents the mean drift velocity of the existing electrons neglecting cre- 
ation or loss of electrons. The bulk coefficients describe the averaged change of the electron 
bulk as a whole. For example, the bulk drift velocity is the displacement of the mean position 
of the electron swarm. The mean displacement of the swarm is not only the result of the elec- 
tron drift, but also of the spatially varying mean electron energies, ionization and attachment 
rates throughout the swarm [60, 9]. Bulk transport coefficients are experimentally measurable 
transport quantities, while the measurement of the flux transport coefficients is rather difficult for 
swarms at high electric field at normal pressure. For more dis cussion of flux and bulk coefficients, 
we refer to Section 2.2. 1 and to [83, 84, 85, 9]. 

In Fig. 4, the transport parameters: electron mobility fi, electron diffusion tensor D, and the 
ionization rate a-, together with the attachment rate a art of electron ensembles are presented for 
both bulk coefficients (marked with "x") and flux coefficients (marked with "o"). In Fig. A.16.e 
and f, the mean energies e and the nonlocal parameter k\ are also presented. For simplicity of 
notation, we fix standard temperature To and vary the pressure p. Quantities really scale with gas 
density no = p/KT. 
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Figure 4: Shown are the transport coefficients, ionization and attachment rate, and the mean energy of electrons in air as 
a function of the reduced electric field at room temperature. The presented coefficients are the bulk coefficients (marked 
with "x") and the flux coefficients (marked with "o") calculated from the particle swarm experiments. Their empirical 
fits are indicated with dashed lines. 
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For the particle swarm generated coefficients presented in Fig. 4, we generated empirical fits. 
They are indicated by dashed lines in Fig. 4. The fitting functions for the bulk coefficients are 



H(E) /m 2 V _1 s- 
ai(E) /m" 
a a „(E) /m" 
D T (E) /mV 
D L (E) /mV 



= exp [-3.18 - 2.65 x 1(T 2 ■ InE + 3A4/E - (4.99/F) 2 ] 
= exp [1.14 x 10 + 3.73 x 10 _1 ■ \n£ - 1.88 x 10 2 /E] 
= exp [l.63 x 10 - 1.95 ■ In E - 8.32 x 10/E] 



5.43 x 10~ 2 + 1.24 xlO -3 



6.01 xl0" 6 -£ 2 + 1.15 xlO" 8 



1 .42 x 10~ 2 + 1 .84 x 10" 3 ■ E - 7A6 x 10~ 6 ■ E 2 + 1 .41 x 10~ 8 



and for the flux coefficients: 



a*(E) /m" 
a* at ,(E) /m" 
D* T {E) /mV 
DUE) /mV 



exp [-2.30- 2.39 x 10" 1 -ln£-6.57x 10 _1 /5 + (6.71 x 10~ l /E) 2 ] 
exp [l. 04 x 10 + 6.01 x 10" 1 -ln£- 1.86 x 10 2 /£] 
exp [l.49 x 10 - 1.63 ■ In E - 7.30 x 10/&] 



5.47 x 10~ 2 + 1.19 x 10" 3 ■ E - 5.06 x 10" 6 ■ E 2 + 1.02 x 10~ 8 



8.87 x 10" 3 +2.26x 10" 3 -E 



8.30 x 10~ 6 -E 2 + 1. 



x 10" 



k x {E)jm = 



6.19 x 10" 7 + 1.89 x 10 _2 /(£ + 4.73 x 10) 2 . 



(19) 



where E = E/(kV cm -1 bar -1 ) is the electric field in dimensionless units, and aj a „ are the ioniza- 
tion and attachment rate. These fitted coefficients will be used in the fluid model to reach optimal 
agreement between particle and fluid model in the hybrid computation. In Appendix Appendix 
A, we also compared our calculation with the Bolsig+ [71, 32] generated transport coefficients 
and reaction rates for air, nitrogen, oxygen, and argon. 

We refer to [4] for the simulation results of the particle model. 

2.2. Fluid model 

2.2.1. The extended fluid model with flux coefficients 

A fluid model for streamers consists of continuity equations for electron and ion densities 
n e , p coupled to the electric field E 



dn e . 
— — + V ■ j e 

dt 



= S, 



drip 



= S, 



Je 

V-E 



-H(E)En e - D(E) ■ V« e 
q (n p - n e ) 



=: -J -J 



(20) 

(21) 
(22) 
(23) 



} e is the electron flux; it consists of an advection part -j fl = - / u(£)E« e and a diffusion part 
-j d = -D(E) • V« f . The source term S accounts for the ionization or attachment reactions that 
change the density of electrons or ions. The ions are approximated as immobile on the time scale 
of electron motion, therefore the total space charge created by different types of positive and 
negative ions can be summarized into one ion density n p . q is the elementary charge, and eo is 
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the permittivity of the gas. The electric field is typically calculated in electrostatic approximation 
E = -V(f>. 

The source term for impact ionization and attachment is traditionally written in local form as 



where the Townsend coefficient a(E) (E — |E|) depends on gas density, collision cross sections 
and the local electric field. It is also tempting to simply insert transport and reaction coefficients 
from electron swarm experiments, i.e., the so-called bulk coefficients from section 2.1.3, into 
the fluid equations. However, in [9] we showed that in high electric fields (above 75 kV/cm in 
nitrogen at standard temperature and pressure) this classical fluid model does not approximate 
the MC particle model well. But when particle and fluid model are inconsistent, discontinuities 
are building up in time at the interface between fluid and particle model when simulating a planar 
streamer ionization front [9]. Therefore, the flux coefficients of section 2.1.3 have to be used, 
and the fluid model has to be extended by writing the source term as 



(Instead of extending the source term, one can also let the source term depend on the mean elec- 
tron energy and introduce another equation for the electron energy, as discussed in the appendix 
of [9]). Although only the source term differs between the classic and the extended fluid model 
and the models describe the same experiments, all transport and reaction coefficients differ be- 
tween the two models. 

There are two reasons to use the extended fluid model with flux coefficients in our hybrid 
model, one is based on the macroscopic simulation results, the other on the microscopic under- 
standing of the electron fluxes. 

First, the extended fluid model approximates the particle model better than the classic fluid 
model. The particle density profiles generated by these three models were compared in [9] for 
electron swarms as well as for planar ionization fronts. The electron and ion density profiles are 
almost identical in the swarm simulations for all three models, but for the classical fluid model 
this is due to a cancellation of terms, as the electrons move too fast, but the ionization rate at 
the tip of the swarm is too low. For the planar front, the electron and ion densities behind the 
ionization front agree well between particle model and extended fluid model, but they are too 
low in the classical fluid model. 

Second, the extended fluid model uses the flux coefficients, that characterize the average 
electron motion. A consistent definition of electron fluxes is important when coupling particle 
and fluid model in space. The definitions of 1) the number of electrons crossing the model 
boundaries in the particle model and of 2) the electron density flux rate at the model boundaries 
in the fluid model have to be consistent. Since the flux coefficients represent the real motion and 
reaction rate of electrons, while the bulk coefficients also depend on the macroscopic electron- 
density-profile, a consistent coupling of particle and fluid model can only be achieved with flux 
coefficients. 

From now on, the term "fluid model" will refer to the extended fluid model with flux coeffi- 
cients. 

2.2.2. Numerical discretization of densities and fluxes 

In our implementation of the fluid model, the electron and ion densities are discretized in 
space with a finite volume method based on the mass balances for all cells, and the electron 



S = \n e ju(E) E| a(E), 



(24) 




(25) 
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Figure 5: Evaluation of the electron flux on the boundary between two cells. 



density is updated in time using a third order upwind-biased advection scheme combined with a 
two-stage Runge-Kutta method. The numerical algorithm and the spatial discretizations of the 
continuity equations in a 2D radially symmetric system have been discussed in [3, 86]. 

We now focus on the numerical discretization of the electron fluxes in the fluid model in 3D. 

The flux of electrons is given in Eq. (22), where the j" = jj(E)En e and j d = D(E) ■ Vn e denote 
the advective and diffusive electron fluxes through the cell boundaries, and D is a tensor in the 
form of 



where I is the identity matrix. 

Electron diffusion in electrostatic fields in gas discharges is frequently approximated as 
isotropic in fluid simulations [37, 42, 39, 87, 88, 89], while it, of course, is anisotropic [90, 91]. 
For example, in the flux diffusion coefficients presented in Fig. 4, the longitudinal diffusion rate 
Dl and transversal diffusion rate Dj relative to the direction of the electric field differ due to the 
anisotropic collision processes in the particle model: Dl < Dj at low fields E < 50 kV/cm, and 
Di > Dj when the field strength is above 50 kV/cm. The discretization of the flux terms requires 
extra care when the diffusion tensor is anisotropic [92]. 

Note that the electron and ion densities calculated at cell centers on a uniform grid, can 
also be viewed as averages over the cell. The electric potential (p and the field strength E are 
taken also in the cell centers, where E determines the electron and ion growth in the cell. The 
electric field components E = (E x , E y , E z ) are taken on the cell vertices, where they determine 
the mass fluxes. The flux coefficients /i and D, and ionization rate a, and attachment rate a a „ 
with anisotropic scattering are from Section. 2.1.3. 

We first consider the 2D example to explain how the electron flux is numerically discretized. 
As shown in Fig. 5, we have a uniform grid with the quantities / evaluated in the centers of the 
grid cells (/ can be n e , <p or E). Consider the flux on the boundary between cell (i, j) and (i+l, j). 
The advective flux /" . can be written as 




(26) 



(27) 
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where n T = (1,0), E^, = (E x ^ i+ 1^, £,,,(,•+ given by Eq. (29), and n ei +ij is approximated 
by a third-order upwind-biased scheme. This gives mass conservation and monotone solutions 
without introducing too much numerical diffusion [86]. The Koren limiter function is used here. 
Denote E L = max(-E, 0) and E R = min(— £, 0) to distinguish the upwind direction for the field 
components, the advective flux j" at cell face (x i+ i ,yj, Zk) then can be rewritten as: 



'x,(f+i,j) 



p(E)E x [n eXiJ) + tff(p UJ )(n e xi+i,j) - ««,(y))] 
+fi(E)E R [n eXi+lJ) + i//(l / p i+ ij)(n eXiJ) - n e ,(i+ij))] 



(28) 



where pij = 1J> , ^ is the limiter function iff(ff) = max (o, max (l, | + f ,#)), and here E 

and E x are taken at cell face {x i+ i, yj). 

To obtain the electron diffusive flux at the cell boundary (x i+ i ,yj), we need the electric field 
and the electron density gradients at the boundary. The electric field is taken as 



-x,(i+i,J) 



1 

AX 
1 

2 



The diffusive flux y* d is calculated as: 



(0(U) - 0a+i,/>) 
2^(*<w-i)-*ft*i)) + 2z;;< 



'(/+1J+1) 



D L (E)- 



EE' 



Dt(E) I 



EE' 



V;; 



(29) 



(30) 



'+5 J 



where the field strength is taken at the cell boundaries and the density gradient V« f at cell face 
(i + i, /) is defined in the same way as the electric field, 



<9n f 

dx (i+ij) 




dn e 


1 


1 


(i+i J) 


2 


2Ay 



( n e,(.iJ) _ n e,(i+l,y)) 



1 



-( n e,(i,j-\) - «e,(ij+l)) + 2^( n e.(»+lJ-l) ~ n e,(i+lJ+l)) 



The flux in the x-direction in 2D therefore can be written as 

je, x = - ;1 = n r ■ (-KE)En e ~ D(E) ■ Vn e ). 



(31) 



(32) 



In the y-direction, Eq. (32) applies as well with x exchanged by y and with n r =(0 1). 

The electron flux on a cell face can be written in the same way for 3D as Eq. (32) for 2D, 
where D is a tensor defined in Eq. (26), n is a vector normal to the cell face, E r = (E x , E y , E z ), 
and (V« e ) r = (dn e jdx, dn e /dy, dn e /dz) are taken at the cell face. For example, the diffusive flux 
is calculated with second-order central differences as 



1 



>x,(i+\,j,k) £2 



dn 



dn e 



= -r\ D L (E) - D T (E) \E X \E x -^+E y -l+E z -l\ + Dj(E)—- 



dx 



dy 



. dn e 



dx 



(33) 



where E, E, and V« e are taken at the cell face (x i+ i, yj, Zk)- Here E at the cell face is defined in 
the same way as in Eq. (29) and V« e is defined as in Eq. (31). 
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Figure 6: The computation time of the Fishpack 3D Poisson solver on m X m X m grids, as a function of m. The test is 
carried out on a desktop computer with AMD Athlon 1.6GHz CPU and 1G memory. 

2.2.3. The Poisson equation 

In the particle as well as in the fluid model, the electric field has to be calculated at each time 
step, therefore this must be done in a very efficient way. A fast 3D Poisson solver fishpack has 
been chosen for this. It uses second-order finite differences. The fishpack is a direct solver, based 
on cyclic reduction with Fast Fourier Transform (FFT) in the third dimension. Although it is not 
versatile (it can only solve Poisson and Helmholtz problems on rectangular grids), it is very fast. 

The computing times are evaluated for a test problem with a Gaussian source term on various 
mxmxm grids. In Fig 6, we plot the computing time as a function of m, when double precision 
is used. Since the solution is direct, the actual shape of the solution does not influence these 
computing times. Another advantage of fishpack is that it needs, apart from arrays to store the 
solution, almost no additional memory. For details and additional tests with this solver we refer 
to [93, 94, 95]. 

2.2.4. Results of a 3D fluid simulation 

Here we present 3D fluid simulation results for a negative streamer in air, propagating in a 
background field of -100 kV/cm. Initially 500 electrons and ions are each distributed within the 
same Gaussian distribution around a spot near the cathode. The simulation was carried out on a 
grid of 256 x 256 X 512 points with Ax - Ay - Az - 2.3 fim, on a system with x e [-0.29 0.29] 
mm, y e [-0.29 0.29] mm and z e [0 1.17] mm. The time step is At = 3 x 10 -13 s. 

In Fig. 7, the electron density (left column), charge density (middle column), and the electric 
field in z-direction (right column) are shown at 0.48 ns (first row), 0.56 ns (second row), and 0.72 
ns (third row). The simulation starts with an avalanche of electrons present initially. After t= 
0.48 ns, a charge layer is clearly formed and the electric field is altered by the space charge. From 
the 0.48 ns to 0.72 ns, the maximal electron density increases from 1.6 x 10 13 to 7 x 10 14 /cm 3 , 
while the maximal field strength increases from 160 kV/cm to about 290 kV/cm. 
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Figure 7: 3D simulation results of the extended fluid model for a negative streamer developing in air in a background 
field of -100 kV/cm. First row: streamer at t = 0.48 ns, second row: streamer at t = 0.56 ns, third row: streamer at t = 
0.72 ns. The columns show from left to right: electron density, charge density, and electric field E z in the z direction. 
Particle densities and fields are represented on two orthogonal planes that intersect with the 3-dimensional structure. 
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3. 3D hybrid model 



3.1. General coupling procedure 

The simulations begin with a few electron and ion pairs followed only by the particle model. 
As new electrons are generated, the number of particles eventually reaches a given threshold, 
after which the simulation switches from the particle simulation to the hybrid approach. If the 
computation time is not a concern, the threshold number can be set as the maximum allowed by 
a real particle simulation. The threshold in our simulation is normally set to 20 million electrons. 

After the number of electrons reaches the threshold, suitable regions to apply the particle or 
the fluid model are found with given criteria to initiate the hybrid computing. To ensure a rea- 
sonable interaction between the models, buffer regions are built by extending the particle regions 
into the fluid region. In the buffer region, the particle movements and the particle densities are 
followed both by the particle model and by the fluid model. 

To update the particle densities and the electric field from one time step to the next, we first 
follow the movement of electrons in the particle region and the buffer region, and the number of 
electrons crossing the model interface is recorded during this time step. The electron fluxes on 
the interface are used as the boundary condition to update the densities in the fluid region. We 
then map the particles in the particle region to the densities on the fluid grid. The electron and 
ion densities are now known in both particle and fluid region, and the new electric field is then 
calculated. With the updated electric field and the particle densities, the new model interface is 
determined. This procedure is repeated in each time step. 

The coupling procedure in 3D has some similarities to the ID coupling. As for planar fronts, 
i) we use the zero-order mapping (see [9]) around the model interface in the particle model to 
avoid particle leaking, and ii) we use forward Euler time stepping instead of two-stage Runge- 
Kutta method in the fluid model to save computer memory and time. However, additional prob- 
lems appear in a real 3D problem. For example, in a planar front, only one point needs to be 
specified to determine the position of the model interface; in 3D, the model interface is a 2D 
surface and extra care is needed; a strongly fluctuating model interface will create large buffer 
regions and dramatically increase the computation cost. The way to count the electron flux over 
the model interface is also different in 3D. In a planar front, the model interface is a straight 
line and all crossing electrons will contribute to the density flux; in 3D, the shape of the model 
interface is more complicated, and crossing electrons need to be defined carefully. 

In the following sections, we consider a particle simulation that starts with the same initial 
condition as the fluid simulation in Section. 2.2.4, and the simulation is carried out on the same 
grid. The particle simulation runs until t * 0.46 ns, when the number of electrons reaches 2 x 10 7 . 
The particle simulation then switches to the hybrid simulation. The problems and our solutions 
during this transition are described in detail in the following sections. 

3.2. First interface construction: column based splitting 

Aiming at following the high energy electrons in a fully developed streamer, fluid and particle 
model will be applied adaptively in suitable regions in the 3D hybrid model. We would like to 
couple the particle and fluid model in such a way that: z) the hybrid model represents the correct 
physics, ii) high energy particles will be included in the particle region, and Hi) the model is 
computationally as efficient as possible. 

The coupling of fluid and particle model was first realized in a planar front as discussed 
in [61, 9]. Two splitting criteria were studied for the planar front. The model interface can either 
be set at the position where the electron density reaches n e = r\ n c mflA when approaching from 
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the nonionized region, or where the electric field is E — % E max {E + in [9]); here n e ^ nax and E max 
stand for the maximal electron density and electric field, and £ and 77 are real numbers between 
(0, 1), and chosen in such a manner that the calculation is efficient and the error is small. 

A 3D streamer can be decomposed into columns with small transversal area parallel to the 
streamer axis. For each column, the previously derived results for the planar front can be used, 
such as the best position of the model interface for a given field ahead of the ionization front. 
When this is implemented numerically, we take a column of mesh cells in the direction of the 
background electric field, and try to locate the position of the model interface. Note that the 
cell-columns near the streamer axis cross through a rather planar ionization front, while the cell- 
columns at the streamer edges are further from this situation. 

3.2.1. Density criteria 

We first try to determine the model interface using the density criterion. But the obtained 
model interfaces have strong fluctuations, especially along the side of the streamer. In Fig. 8(a) 
and 8(b), we show the position of the model interface when it is set at n e - n ejnax (left) or 
n e = 0.7 n e , max (middle). The fluid model is applied at regions above the model interface and the 
fluid model is applied beneath. The model interface is set at n e = n ejmax or n e = 0.7 n emax only 
when n ejnax > c, where c is the density value for one electron per cell. It means that we split a 
cell-column into particle and fluid region only if it contains electrons; in the large area without 
electrons, the particle model is applied. When the model interfaces are placed at n e = n e%max , 
the fluctuations are almost everywhere. When the model interfaces are placed at n e = 0.7 n eylmlx , 
the model interface is a rather smooth surface in the center where n emax is relatively large, and it 
fluctuates strongly where n efnax is small. 

The density fluctuation is the first problem when we apply the planar front coupling to the 
cell-columns in 3D. The electron density fluctuations in a planar front can be suppressed by ex- 
panding the planar front in the transversal direction, which increases the number of electrons and 
smoothes the density profile. This method can be applied for the planar front since the parti- 
cle distributions are uniform in the transversal direction and therefore charged particle densities 
remain unchanged. But this numerical trick can not be applied for the cell-column in 3D cou- 
pling, while the limited number of electrons or ions within the cell-column gives rise to strong 
fluctuations in density profiles along the column. 

To reduce the computational cost and save memory for high energy electrons at the streamer 
head, we would like to have the model interfaces of all cell-columns to lie within a smooth 
surface. A fluctuating model interface requires extra buffer regions and therefore more electrons 
to be followed by the particle model, and it creates a large volume of the model interface in the 
transversal directions where the motions of all electrons close to the lateral model interface have 
to be traced. It is therefore not a good idea to relate the position of the model interface to the 
maximum density within a column. 

When the interface is placed at smaller densities, e.g., at cells at position n e = 0.7 n e ^ max , it is 
rather smooth over z, at least in the center of the streamer where n e<max is large enough. For planar 
fronts in fields between 50 and 200 kV/cm, the particle densities from a hybrid computation with 
a model interface at n e = 0.7 n ejnax deviates from the pure particle simulation results by not more 
than 4% [9]. 

To obtain a smooth model interface at the edge of streamer, one possible solution is to find 
a fitting formula for the interface position near the streamer center and to extrapolate it to the 
side to calculate the interface position even if no electron exists in that cell-column. But besides 
concerns on the asymmetry of the density over the radius and on the difficulties of fitting, the 
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(a) n e = Uejnax 



(b) n e = 0.7 n, 



(c) E = 0.84 E,„ 



Figure 8: The position of model interfaces for a streamer at t»0.46. The model interface is at n e = n ejnax (left), 
n e = 0.7 n ema (middle) and E = 0.84 E max (right) along each cell-column. Far outside of the streamer where no free 
electrons and ions exist, the particle model is applied and the model interface is set at z = 0, When the criterion is 
n e = n ejnax , fluctuations are everywhere. When n e = 0.7 n emax the interfaces are rather smooth tn the center of the 
streamer, the fluctuations appear at the side of streamer. When E = 0.84 E max , the model interface forms a smooth 
surface. 

interface position at the streamer side would not depend on local properties in this procedure, but 
on the fit to the streamer center. 

3.2.2. Field criteria 

In a planar front, there is a clear correlation between the local electric field and the electron 
density, and the electric field is always smoother than the electron densities. Therefore in the 
planar hybrid calculation, also an interface criterion depending on the field level E = i; E max 
was tested. The relation between particle density and electric field in 3D is not as clear as in 
ID. Therefore a field dependent criterion is tried with several values for the particular value 
of £ = 0.84 the positions are given in the Fig. 8(c). The field criterion £ = 0.84 for the model 
interface agrees very well with the density criterion rj = 0.7 in the center of the streamer, while 
model interfaces are much smoother at the side of streamer. 

One problem of the field criterion is that at the side of the streamer where the field varies less 
than in the center, it becomes E < 0.84 E max everywhere along the cell-columns. Because there 
are only few electrons which occasionally fly out of the channel in the sidewards direction, we 
leave them to the fluid model. Therefore from the point where the level E = 0.84 E max ceases to 
exist, the model interface is extended horizontally several cells outwards to include all particles 
at the streamer side into the fluid region. 

At the streamer side, we also add one or two extra cells for the density diffusion in the fluid 
calculation. In the particle model, the electrons are discrete and the electron density is always 
Nxc, where c is the density value for one electron per cell. But in the fluid calculation, diffusion 
can fill the entire fluid region with a small n e > 0. To avoid a continuous expansion of the fluid 
region to the side, the fluid electron density is set to zero if it drops below c. 

3.2.3. After splitting 

Having described this procedure for locating the model interface, we now present the splitting 
results when a 3D simulation is transferred from the particle calculation to the hybrid calculation. 
In Fig. 9, the electron density (upper left), charge density (upper right), electric field (lower left) 
and electron mean energy (lower right) at t« 0.46 are shown with the model interface and the 
buffer region marked with red lines (blue lines for the electron mean energy). Below the model 
interface, the fluid model is applied; above the model interface, the particle model is applied. 
With the cell-column based approach, the particle model focuses on the electrons at the streamer 
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Figure 9: Electron density (upper left), charge density (upper right), electric field (lower left) and the electron temperature 
(lower right) at tw 0.46 are presented with the model interface and buffer region marked. 

head where the local electric field is most enhanced. And it leaves the large remaining region 
with high electron densities and low electric field to the fluid model. 

The first results of 3D hybrid calculation results based on this approach were reported in the 
letter [11] for a negative streamer in N2. The particle region in this approach focuses on the 
electrons at the streamer head and leaves the region behind it to the fluid model. However in the 
situation where i) more than one streamer propagates in the simulation domain and the streamers 
are not in parallel, ii) a runaway electron runs out of the planar front and creates an avalanche 
ahead of the streamer, or iii) a double headed streamer emerges, this cell-column based approach 
is unable to deal with them. To deal with those situations, a more flexible splitting algorithm is 
needed. 

3.3. Second interface construction: full 3D splitting 

We have developed another splitting method that is simultaneously based on the electric field 
and on the electron density. In contrast to the previous splitting method that was based on local 
column-based quantities, we now develop a global splitting criterion that operates on the whole 
3D domain. 

The particle regions are now taken as the regions where the electric fields are higher than a 
global field threshold and where electron densities are lower than a global density threshold, and 
the fluid regions covers the rest. Most of the energetic electrons, and electrons with the potential 
to gain high energies are in the region with the strongest electric field. 

In Fig. 10, we show the model interfaces (marked with red lines) determined only by a field 
threshold. The presented streamer is from the same simulation at same time t« 0.46 as shown 
in Section 3.2. The particle model is applied if the electric field is strong enough, E > ifE/,, in 
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Figure 10: The model interface based on the electric field only. Electron density (left), charge density (middle), and 
electric field (right) at t« 0.46 are presented with the model interface marked. The particle model is applied at the region 
where the E > f £7, with £ = 0.7 to 1.0, and the fluid model is applied to the rest. 
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which Ei, is the background field 100 kV/cm and ^ is a free parameter. Note that we are not 
searching a specific model interface anymore, we now search for regions suitable for the particle 
model. Fig. 10 shows how the particle and the fluid region change when ^ changes from 0.7 to 
1 .0. The fluid region has the smallest volume with £ = 0.7 and the largest volume with f = 1.0. 
From £ = 0.7 to 1.0, the number of electrons in the particle region decreases from 9.6 x 10 6 to 
2.5 x 10 6 . 

Thresholds with £ < 0.9 are not appropriate to determine the model interface, not only be- 
cause there would be too many electrons within the particle region, but also due to the fact 
that electrons at the side of the streamer can cause a considerable increase of the computational 
cost when they cross the model interface. In both cases £ = 0.9 and 1 .0, electrons at the side 
of streamer are covered by the fluid regions. But £ = 1.0 leaves all the lateral region to the 
fluid model which makes the model less flexible and computationally more expensive. The fluid 
model is computationally more efficient than the particle model only in mesh cells that contain a 
considerable amount of electrons. Letting the fluid model cover the large open side regions will 
increase the computational cost due to the non-zero, but sometimes unphysically small densities 
at side. 

We tested a range of interface positions from £ = 0.9 to £ = 1.0 and the field criterion 
£, = 0.94 was chosen for the hybrid simulation. With this threshold, the sides of the streamer are 
fully covered by the fluid model during the period of interest. Together with the field criterion, 
a density criterion is added. The fluid model is applied if the electron density is high enough, 
n e > r\ n ejnax , where n e>max is the maximal electron density in the whole simulation domain. In 
Fig. 1 1, we present the electron density, the electric field and the local mean energy of electrons 
with the model interface indicated with a curved line. The particle model is applied in the region 
where the electric field E > fEj, and electron density n e < rj n e max , and the fluid model is applied 
in the region where the electric field E < £Eb or electron density n e > rj n ejnax , where f = 0.94 
and rj = 0.7 (upper panel) and 0.9 (lower panel). 

As shown in Fig. 10, since the side of the streamer has already been allocated to the fluid 
region, differences only appear at the streamer head due to the choice of rj in the density criterion. 
When the density threshold changes from 77 = 0.9 to a small value such as 77 = 0.7, the particle 
region shrinks and the fluid region expands towards the field enhancement region. 

In contrast to the column based splitting, which can only follow one streamer head, the full 
3D splitting approach can apply the particle model in separate regions. For example, it can follow 
both the negative and positive ionization fronts when one simulates a double headed streamer, or 
a second streamer head if another streamer forms. 

Hybrid simulations for a double-headed streamer in air with photo-ionization have been car- 
ried out in an overvolted gap. The photo-ionization model follows Zheleznyak's model [96], 
and the implementation for the particle model is the same as in [26, 97]. Although the devel- 
opment of the double headed streamers, and multi-streamers are not the topic of this paper, two 
time steps from such simulations are presented in Fig. 12 to show the flexibility of the full 3D 
splitting approach. 

Fig. 12 shows a double-headed streamer which is developing in the middle while several 
small avalanches appear around the streamer. On the left we plot the equal-density surface of 
the electron bulk with electron density n e > 10 13 /cm 3 . On the right we plot the electric field on 
two orthogonal planes intersecting the 3D structure, with the model interface and buffer regions 
marked with red lines. In the left plot of the upper panel, the larger electron bulk is the main 
streamer, while some other small avalanches develop around it. The small avalanches are caused 
by the electrons created by photo-ionization; i.e., by photons emitted from the streamer head, 
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77=0.7 



77=0.9 




Figure 1 1 : The location of the model interface is simultaneously based on a field and a density threshold. Electron 
density (left), electric field (middle), and electron mean energy (right) at t» 0.46 are presented with the model interface 
marked. The particle model is applied where the electric field E > 0.94£j and the electron density n e <r\ n ejnax with 
r\ = 0.7 (upper panel) and 0.9 (lower panel), and the fluid model is applied in the remaining region. 



t = 0.51 ns 



t = 0.6 ns 




Figure 12: A double headed streamer is developing in air while photo-ionization is included. A surface of equal electron 
density (left) and the electric field in color coding (right) are presented; the model interface with the buffer region are 
marked. 
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and several of them have already created a relatively large electron density reaching the density 
threshold. Separate fluid regions therefore are created in this area. It happens that the xz plane 
where the electric field is plotted crosses one of these avalanches. A cross section of the separate 
model interface for this avalanche is therefore also included in the electric field plot. 

In the lower panel, we plot the streamers 0.09 ns later. The separate fluid region for the 
small avalanche has merged with the main fluid region of the main streamer channel, while a 
new separate fluid region has formed. Note that not only the upper streamer head is followed 
by the particle model, but that the particle model also follows the field enhanced region of the 
downward propagating streamer head. 

3.4. Structure of the buffer region 

The interaction of two models is realized through the well know "buffer region" technique 
which has been employed in hybrid computations for air flow [98, 99, 100, 101, 102, 103, 104], 
liquid flow [105, 106, 107], and also in small scale solid systems [108, 109, 110]. A similar 
method that is more suitable for streamer simulations has been implemented and tested for planar 
fronts [61,9]. Here our implementation in the 3D hybrid model is discussed. 

The electrons that cross from the particle region over the model interface will be recorded, 
since they contribute to the density in the fluid region. We have explained in [9] that in a planar 
negative front, since the electrons move on average slower than the model interface, we only 
need to remove the electrons from the particle list if they fall into the fluid region, but we do not 
need to artificially create new electrons in the buffer region to stabilize the electron flux at the 
model interface. The same occurs in the (negative) head of a 3D streamer: electrons on average 
move slower than the model interface in the propagation direction, and we don't need to add 
new electrons artificially. The side of the streamer is normally attributed to the fluid region in 
our 3D hybrid model; as no electrons are followed by the particle model here, the interface does 
not need to be constructed here. At the tail of the streamer, there is no electron flow without 
photoionization or background ionization, and the model interfaces are more or less stationary. 
If photo-ionization is added, at the positive head of the streamer, the electrons propagate into the 
streamer from the particle into the fluid region. Therefore no new electrons need to be created 
artificially at any part of the model interface. 

Here we set the length of the buffer region as 3 cells in the z-direction and 2 cells in x- and 
y-direction with the cell length Ax = Ay = Az — 2.3 yum, since the maximum electric field in z- 
direction is normally much higher than in the x- or y- direction (see Table I in [9]). 

It is important in 3D that for any cell face, cell edge and cell corner shared by the particle 
region and the fluid region, a buffer region separates them. A direct contact of particle and fluid 
model without a buffer region can cause electron leaking, which creates loss of mass and charge. 

The so called "corner problem" [100] is a technical but important issue when calculating the 
electron flux in the hybrid model. When an electron passes from the particle region into the fluid 
region or vice versa, its contribution to the flux on the cell face is recorded if the cell is in the 
fluid region, since only the density update in the fluid region needs the electron flux. The cell 
face is not necessarily the model interface between particle and fluid region. A 2D cartoon is 
shown in Fig. 13 to illustrate the problem. The particle region is in the upper right part and the 
fluid region is in the lower left corner, the particle region extends into the fluid region by 2 cells 
in all directions and creates a buffer region. As illustrated in Fig. 13, when an electron flies from 
cell "al" to "bl", it contributes to the flux on the model interface at cell face "al<-> bl". A more 
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Figure 13: Schematic illustration of the corner problem between the particle and fluid cells around the model interface. 
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e (eV) 

Figure 14: Electron energy distribution function (EEDF) of all electrons in the streamer at t»0.46 ns (solid line) and the 
part of electrons in the particle and buffer region (dotted line). The EEDF is not normalized and Ae = 0.2 eV. The part of 
electrons which removed from the particle list are mainly the low energy electrons. 



complicated case occurs when an electron flies from cell "b2" to "cl", it contributes not only to 
the model interface at cell face "b2«-> c2", but also to the cell face "cl«-> c2". On the other hand, 
when the electron flies from cell "c2" to "b3", it only contributes to the model interface at cell 
face "b2 <-> c2". Finally, consider that the electron flies from "b3" to "c4"; it flies over two model 
interfaces, but in the model it contributes to neither of them since there is no mass change in the 
fluid region. 

When the simulation reaches a total number of 20 million electrons, at t^0.46 ns, we split the 
simulation domain into a particle and a fluid region. If we use the column based splitting with 
model interface at E — 0.84 £,„, 12 million electrons which are neither in the particle nor in the 
buffer region are removed from the particle list and transferred into particle densities in the fluid 
region, while 5 million electrons remain in particle region and 3 million electrons remain in the 
buffer region. If we define the model interface using the full 3D splitting with E = 0.94 E\, and 
n e = 0.7« e mnv , 15.5 million electrons in the fluid region are removed from the particle list, while 
2.5 million electrons remain in particle region and 2 million electrons remain in the buffer region. 



26 



t=0.48 



t=0.56 



t=0.72 




electron density (1 /cm 3 ) Electric field (kV/cm) High energy electrons (eV) 

Figure 15: 3D simulation results of the 3D hybrid model for a negative streamer developing in a background field of -100 
kV/cm. First row: the start of hybrid computation at t = 0.48 ns, second row: streamer at t = 0.56 ns, third row: streamer 
at t = 0.72 ns. The columns show from left to right: electron density, charge density, and electric field strength. Particle 
densities and fields are represented on two orthogonal planes that intersect with the 3-dimensional structure. 

Both splittings leave most of the ionization front to the particle model. It is also remarkable that 
the majority of the high energy electrons remain in the particle list which results in a good model 
for the study of runaway electrons. In Fig. 14, we show the electron energy distribution function 
(EEDF) of electrons from the whole streamer and from only the particle and buffer region. The 
EEDF is not normalized so one can clearly see in which energy range the electrons are removed 
from the particle list. 

4. Simulation results 

Having introduced the new coupling scheme with the numerical details, we now present our 
hybrid simulation results for streamers in air without photo-ionization. 

The initial conditions and the configuration are the same as the fluid simulation in Sec- 
tion 2.2.4. The model interface is determined by the full 3D splitting approach with the field 
criterion Eg = Q.94E/, and density criterion n ejj = 0.1n e>max . In Fig. 15, we show a negative 
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streamer followed by the hybrid model at different times. At time 0.48 ns (upper panel), after the 
hybrid calculation started 0.02 ns, there are 2.8 x 10 7 electrons in total, where 6.0 x 10 6 electrons 
are followed by the particle model in the particle and buffer region, and the maximal electric field 
at that time is around 150 kV/cm. 

At time 0.56 ns (middle panel), the electric field reach 180 kV/cm, there are 7.5xl0 7 electrons 
in total, in which 8 x 10 6 electrons followed by the particle model. 

In the bottom panel of Fig. 15, we show the streamer at t=0.72 ns. The maximal electric field 
at that time is around 300 kV/cm. There are 4.5 x 10 8 electrons in total, and 2.0 x 10 7 electrons 
are followed by the particle model. 

At time 0.56 ns, electrons with energies above 200 eV start to appear at the streamer head. 
But electrons with such energies can not hold their energy long, they lose their energies very 
quickly and disappear in the next time step. As the streamer propagates, the maximal electric 
field keeps increasing, there are more electrons with energies above 200 eV generated and the 
maximal electron energy also increases. At t= 0.72 ns, the maximal electron energy is around 1 
keV. 

5. Conclusion 

5.1. Numerical technique 

The new splitting method, the full 3D splitting, has overall advantages over old column- 
based splitting. It adaptively applies the particle model in the most field-enhanced region and the 
fluid model in regions where the density is high enough. It is flexible in complicated situations 
where spatial coupling is required. Moreover, it is easier to add new features such as local grid 
refinement in the full 3D splitting than in the column based splitting. 

The accurate definition and the evaluation of the electron transport coefficients and reaction 
rates from the particle model for the fluid model are essential for the subsequent successful 
coupling of the models. This is the reason why a large part of this paper and also our previous 
papers [60, 9] are devoted to the correct calculation of these coefficients. These coefficients are 
required for the successful implementation of the hybrid model. 

Compared to simulations with the pure fluid or pure particle model, the 3D hybrid model 
combines the advantages of both models. It describes the dynamics of a streamer channel in a 
very efficient way while being able to follow the movement of each single electron in the most 
active region of the streamer head. The 3D hybrid model is therefore a powerful tool to study the 
kinetics of electrons in the important regions of streamers . 

5.2. Physical predictions 

The 3D hybrid streamer model is a major step forward for several physical problems. First, it 
correctly traces the high electron energies in space and therefore can supply correct densities for 
the many different excited states of molecules; these distributions determine the further chem- 
ical reactions and end products at later times. Second, the model traces both electron density 
fluctuations and electron energy fluctuations, and therefore it can correctly trace the influence of 
fluctuations, e.g., on streamer branching. These two aspects have not been evaluated yet from 
the simulations. Third, as shown in [11], it can predict electron run-away from streamers. We 
now elaborate on this aspect. 

Electrons with energies above 200 eV are shown in Fig. 15 and similarly in [11]; they appear 
in the regions with high local field enhancement, and some of them are accelerated up to 2.5 keV 
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before they disappear into the anode. The energy of 200 eV is the threshold value for electron 
runaway in air. As shown in Fig. 2, above 200 eV, the electron collision frequency decreases, 
and the electrons are accelerated further more easily. 

The development of run-away electrons has been studied in a range of constant uniform 
electric fields with particle swarm experiments [76, 24, 111, 70, 14, 15], in agreement with [112, 
10], where it was shown that thermal electrons run away when electric field exceeds a critical 
strength of E« 260 kV/cm at standard temperature and pressure, and the runaway rate increases 
as the field increases. This is quite reasonable if one considers that in such a field an electron 
gains 200 eV within less than 8pm, which is the range of the mean free ionization length. 

However, much lower background fields together with the field enhancement at the streamer 
tip are sufficient to let electrons run away. In our simulation, the background field is 100 kV/cm 
and the streamers are quite short (as the simulations do not have adaptive grid refinement yet). 
The field in these simulations is eventually enhanced by a factor of 3, but even when it just 
exceeds 180 kV/cm locally, electron run-away sets in. Chanrion et al. even observed electrons 
with 1 keV when the locally enhanced field reaches 160 kV/cm N/No for a negative streamer in 
air at 70 km altitude [12] (where A^o is air density at ground level and N at 70 km altitude). The 
discrepancy between our and Chanrion's data is probably due to the fact that the total number of 
electrons in a streamer scales as 1 /N [2], and therefore streamers at high altitudes contain more 
electrons, and therefore the rate of run-away electrons is larger. But how can the different run- 
away thresholds of 260 kV/cm in homogeneous fields and of 160 or 180 kV/cm in the streamer 
head be explained? 

The answer is two-fold. First, 260 kV/cm is the threshold where the majority of electrons 
experience more acceleration than friction and run away, while single electrons in the high energy 
tail of the distribution can run away earlier. Second, the electrons in the streamer head are not in 
equilibrium to the local electric field. Even in a constant field, the electrons in the tip of a swarm 
on average have higher energies than at the back end [60, 9]. It is these high energy electrons at 
the tip of the ionization front that run into the highly enhanced electric field at the streamer tip 
and are accelerated further. 

But if the electrons run away, they run into the region with lower field ahead of the streamer. 
For example, electrons with an energy of 200 eV run ~6 times faster than the streamer when its 
maximally enhanced electric field is 250 kV/cm. Depending on the spatial profile of the electric 
field and on the energy of the electrons, they will either predominantly loose their energy ahead 
of the streamer or accelerate further (if no anode is in their way). 

For fully understanding the hard radiation from sparks [20, 21, 22] and corona streamer 
discharges [23], two ingredients are necessary: the acceleration of thermal electrons to energies 
above 200 eV in streamers or next to pointed electrodes; and the profile of the electric field 
ahead of a streamer corona or next to another strongly curved electrode that supports the further 
acceleration of the energetic electrons. 

5.3. Outlook 

The size of the system is kept small to obtain sufficient numerical accuracy, while accurate 
solutions are not required everywhere. In the future local grid refinement will be incorporated 
in the hybrid model. It will allow us to study streamer propagation in a large system, and the 
generation rate of the runaway electrons can therefore be obtained even for lower background 
fields. 

We have presented preliminary results for a double headed streamer in air with photo-ionization 
in Section 3.3. Photo-ionization is a three-step process where first a nitrogen molecule is excited 
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by electron impact, then it emits a photon, and the photon ionizes oxygen. In our particle model 
the nitrogen impact excitation is already included. Therefore it will be possible to develop the 
widely-used photo-ionization model further to take the lifetime of the excited species into ac- 
count. It may change the position of the photon source and may consequently influence the 
streamer propagation. 
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Appendix A. Comparing transport and reaction parameters of the fluid model for air, ni- 
trogen, oxygen, and argon 

The electron transport coefficients, attachment and the ionization rate in air have been pre- 
sented in section 2.1, The simulated air is composed by 78.12% nitrogen, 20.946% oxygen , and 
0.934% argon. In this section the generated transport coefficients and reaction rate are compared 
with the bolsig+ package [71]. bolsig+ is a Boltzmann solver to calculate electron transport 
coefficients in gases or gas mixtures. It is based on the two-term Legendre expansion solution of 
the Boltzmann equations [113, 32]. And the bolsig+ calculated transport coefficients have been 
widely used as an input of the fluid simulation for various plasma applications. 

To compare with each other, both our particle model and the bolsig+ package use the same 
cross section data (siglo database [71]) with the same energy splitting mode in ionization colli- 
sions. The energy splitting between the primary and the secondary electrons during an ionization 
has two different modes in bolsig+, "equal sharing" or "primary electron takes all". In both 
bolsig+ and our particle model, we used "equal sharing". 

Both the bulk coefficients and the flux coefficients are calculated from our particle model. 
For each of them, the particle swarm experiments are carried out both with isotropic and with 
anisotropic scattering. Since the momentum transfer cross section is fixed, as discussed in Sec- 
tion 2.1.2, the electron swarms shall be similar in both cases. 

Appendix A.l. Air 

In Fig. A. 16, the transport parameters: electron mobility p, electron diffusion rate D, mean 
energies e, and the ionization a, together with the attachment rate a att , and the nonlocal parameter 
k\ of electron ensembles, are presented. 

They are generated from bolsig+ (solid line), the flux coefficients from our particle simulation 
with isotropic scattering (marked with "□"), flux coefficients with anisotropic scattering (marked 
with "o"), the bulk coefficients with isotropic scattering (marked with "+"), and bulk coefficients 
with anisotropic scattering (marked with "x"). The presented transport parameters from bolsig+ 
have the same definition with the flux transport coefficients. We also made the empirical fittings 
for both flux and bulk coefficients from our particle swarm experiments and they are plotted with 
dashed lines. 

As shown in Fig. A. 16, the coefficients with isotropic or anisotropic scattering have no differ- 
ence for low electric field (E < 50 kV/cm). In the high field range (E > 100 kV/cm), the electron 
swarm simulated with anisotropic scattering has slightly larger mean energies and slightly higher 
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Figure A. 16: Shown are the transport coefficients, ionization and attachment rate, and the mean energy of electrons 
in air. The presented coefficients are the bulk coefficients and the flux coefficients calculated from the particle swarm 
experiments with both isotropic and anisotropic scattering, and with equal energy sharing in the ionization collision. The 
results are compared with results of the Boltzmann solver Bolsig+ [71, 32]. 
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mobilities than that of isotropic scattering; the transversal diffusion rate is larger with isotropic 
scattering and the longitudinal diffusion rate is larger with anisotropic scattering; the ionization 
rates are higher with anisotropic scattering but almost undistinguishable. In general, the trans- 
port coefficients obtained with isotropic scattering are very close to the one with anisotropic 
scattering. The flux coefficients from our particle swarm simulations is undistinguishable to the 
Bolsig+ calculated coefficients for the low field, while small differences remain at strong electric 
fields. 

The fitting functions for the bulk coefficients with the anisotropic scattering are 

H(E) /m 2 V -1 s -1 = exp [-3.02- 6.22 x 1(T 2 ■ ln£ + 2.19/.E - (1.67/£) 2 ] 

<*,(£) /nT 1 = exp [1.13 x 10 + 3.81 xl0 _1 -ln£- 1.85 x 10 2 /£] 

a a „(E) /nT 1 = exp [l. 63 X 10 - 1. 94- In E- 8.24x10/5] 

D T (E) /mV 1 = 5.34xl0~ 2 + 1.35 x 10~ 3 ■ £ - 6.52 x 10~ 6 ■ £ 2 + 1.38 x 10~ 8 ■ £ 3 
D L (E) /mV 1 = 6.83 x 10~ 3 + 2.40 x 10~ 3 ■ E - 1 .00 x 10~ 5 ■ E 2 + 2.23 x 10~ 8 ■ E 3 

and for the flux coefficients: 

lf{E) /n^vV 1 = exp [-2.32 - 2.36 x 10" 1 ■ ln£ - 4.83 x 10 _I /£ + (2.94 x lO" 1 //?) 2 ] 

a*(E) /irT 1 = exp [l.03 x 10 + 6.25 x 10"' \n£- 1.80 x 10 2 /e] 

a* att (E) /nT 1 = exp [l.50 X 10 - 1.65 • \nE - 7.35 X 10/e] 

D* T (E) /mV 1 = 5.43 x 10~ 2 + 1.31 x 10~ 3 ■ £ -6.69 x 10~ 6 ■ £ 2 +6.85 x 10~ 8 ■ £ 3 
£>*(£) /mV 1 = 1.33 xl0~ 2 + 1.95 x 10~ 3 1.17 xlO~ 5 -.E 2 + 2.76 xl0~ 9 -.E 3 
ki(E) /m = 5.58 x 10~ 7 + 2.94 x 10~ 2 /(£ + 8.44x10). (A.l) 

We recall that the fitting functions for air in Section. 2.1.3 are from particle swarm simulations 
with Opal's formula applied to distribute the electron energies between two out-comming elec- 
trons in an ionization, while the functions presented here using an even spliting between two 
out-coming electrons in ionization. 

Although more or less the same swarm and transport parameters can be obtained for the 
different scattering method, the coefficients from the anisotropic scattering are preferred in our 
hybrid model. Not only because anisotropic scattering is closer to physical reality, but also 
because the small variance in the electron mean energy at relatively high fields can make a 
difference for the presence of high energy electrons. The choice of the differential cross sections 
has direct influence on the electron energy distribution function (EEDF). In Fig. A. 17, we show 
the EEDF in a particle simulation with both isotropic (curve marked with "+") and anisotropic 
(curve marked with "*") scattering. The two EEDFs in general agrees with each other very well. 
However, when we zooming into the low energy part 0-15 eV and slightly higher energy part 20- 
45 eV of the EEDF, the comparison clearly shows that there are more low energy electrons with 
isotropic scattering and high energy electrons are easier to produce with anisotropic scattering. 
That is, the scattering method influences the generation rate of the high energy electrons. Since 
one goal of the hybrid model is to study the generation of runaway electrons, the anisotropic 
scattering is used in our particle model, and the parameters from the anisotropic scattering are 
used in our fluid model. 
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Figure A. 17: The electron spectrum in the swarm experiments at 200 kV/cm. The "+" marked curve shows the parti- 
cle simulation results with anisotropic scattering. The "x" marked curve is obtained from the particle simulation with 
isotropic scattering. Two subplots zoom into the 0-15 eV and 20-45 eV of the EEDF. The zoomings show the influences 
of the different scattering methods: more low energy electrons appeal' when isotropic scattering is used, while more high 
energy electrons are generated when the anisotropic scattering is used. 



Appendix A.2. Nitrogen 

The electron transport coefficients and the ionization rates in nitrogen have been reported 
in [60], and the parameter k\ has been reported in [9]. However, the previous particle swarm 
calculation fixed the total elastic cross sections and didn't take into account the corresponding 
changes when different scattering methods are applied. The previously calculated parameters 
with isotropic scattering and anisotropic scattering therefore have noticeable differences, espe- 
cially at electric fields above 150 kV/cm. Now the momentum elastic cross section is fixed, 
therefore we here present a new set of the transport parameters and reaction rates for nitrogen. 

In Fig. A. 18, the transport coefficients p, D, ionization rate a,-, nonlocal parameter ki, and 
the mean energy of electrons e are presented for nitrogen, in which the p, a,, Dj, and e are also 
compared with the solution of bolsig+. As for air, we were also able to find empirical fittings for 
the parameters in nitrogen. 

The fitting functions for the bulk coefficients are 

p(E) /m 2 \- l s- 1 = exp [-3.79 + 8.25 x 10~ 2 ■ ln£ + 8.24/£ - (2.21 x 10/£) 2 ] 

a t (E) /nT 1 = exp [l. 17 x 10 + 3.06 x 10"' -In/? -2.11 x 10 2 /E] 

D T (E) /mV 1 = 4.85 x 10~ 2 + 1.10 x 10~ 3 ■ £ -3.67 x 10~ 6 ■ £ 2 +6.78 x 10~ 9 ■ £ 3 
D L (E) /mV 1 = 7.56 x 10~ 4 + 2.24 x 10" 3 ■ E - 7.73 x 10~ 6 ■ E 2 + 1 .89 x 10 -8 ■ £ 3 
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Figure A. 18: Nitrogen: transport coefficients, ionization rate, the mean energy of electrons, and ki presented as in 
Fig. A. 16 for air. 
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and for the flux coefficients: 



H*(E) /nrV-'s" 1 

a*(E) /m" 1 

D* T (E) /mV 1 
D* L (E) /m 2 s- 1 
h{E) fm 



= exp[-2.58- 1.84 x 10" 1 ■ \nE - 9.60 x IQ~ 1 /E + (4.11/£) 2 ] 



exp[l.05 



x 10 + 5.83 x 10" 1 InE- 2.05 x 10 2 /E 



= 4.79 x 10 - + 1.16 x 10" J ■ E - 4.65 x 10" b ■ E L + 8.12 x 10" 
= 9.08 x 10" 3 + 1.72 x 10" 3 ■ E - 5.90 x 10" 6 ■ E 2 + 1.03 x 10" 
= 4.90 x 10" 7 + 7.58 x 10" 2 /(£ + 1.80 x 10 2 ) 2 . 



(A.2) 



Appendix A. 3. Oxygen 

The electron transport coefficients in oxygen are shown in Fig. A. 19. The electron mobilities 
p., ionization rate a, and attachment rates a att , transversal diffusion Dj, longitudinal diffusion 
Dl, mean energy e, and the coefficient k\ from Eq. (25) are presented in Fig. 19(a), Fig. 19(b), 
Fig. 19(c), Fig. 19(d), Fig. 19(e), Fig. 19(f) respectively. As for the air and nitrogen, we compared 
our coefficients with the calculation results from bolsig+ for mobilities, reaction rate, transver- 
sal diffusion and the mean energies, and a rather good agreement have been obtained between 
the flux coefficients calculated from our particle swarm simulation and the bolsig+ except the 
diffusion rates. 

The difference between transversal diffusion rate Dj with the isotropic and anisotropic scat- 
tering are already remarkable at the very low electric fields . Similar behaviors have been also 
observed in Dj in argon and it is not clear the reasons causing this. Empirical fittings for the 
fluid parameters are given in the following. The fitting functions for the bulk coefficients are 



p(E) /m 2 V"V 
a t {E) ImT 
a att {E) /m" 
D T (E) /mV 
D L (E) /mV 



= exp [-1.75- 3.02 x 10" 1 -ln£-4.34x lO^/E - (8.68 x lO^/F) 2 ] 
= exp[l.ll x 10 + 4.52x10"' -ln£- 1.44x 10 2 /E] 
= exp [l. 26 x 10-9.71 x 10"' \n£- 1.53 x 10/E] 



1.56 x 10" 



1.12 x 10" 4 • E + 1.90 x 10" 6 • E 2 - 4.24 x 10" 9 • E 3 



1 .30 x 10" 1 + 1 .46 x 10" 4 ■ E + 2.99 x 10" 6 • E 2 



8.29 x 10" 9 -E 3 



and for the flux coefficients: 



p\E) /m 2 V _1 s- 

a*{E) /m" 

a att (E) /m" 

D* T (E) /m 2 s~ 
D* L {E) /m 2 s" 
h{E) /m 



= exp [-1.01 -4.89 x 10"' ■ ln£ - 3.16/£ + (3.26/F) 2 ] 

= exp[l.01 x 10 + 6.87 x 10" 1 -ln£- 1.39 x 10 2 /£] 

= exp [l. 19 x 10 -8.10 x 10" 1 - InE - 1.32 x 10/E] 

= 1.52 x 10" 1 + 1.17 x 10" 4 • E — 1.04 X 10" 6 ■ E 2 + 3.35 x 10" 9 • E? 



1.44 x 10" 



1.18 x 10" 4 • E + 1.34 x 10" 5 • E 2 - 3.94 x 10" 8 ■ E 3 



4.17 x 10" 7 + 2.36 x 10" 2 /(£ + 7.07 x 10) 2 . 



(A.3) 
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jure A. 19: Oxygen: transport coefficients, ionization and attachment rates, and the mean energy of electrons presented 
in Fig. A. 16 for air. 
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gure A. 20: Argon: transport coefficients, ionization and attachment rates, and the mean energy of electrons presented 
in Fig. A. 16 for air. 
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Appendix A.4. Argon 



The same parameters and comparisons for argon can also be found in Fig. A. 20, and their 
empirical fittings are listed here. For the bulk coefficients, we have 

H(E) /nrV-V 1 = exp [-3.87 + 7.69 x 1(T 2 ■ In E + 2A2/E - (1.66/F) 2 ] 

a t (E) /m" 1 = exp [9.07 + 7.29 x 1CT 1 ■ In £-3.56 x 10/e] 

D T (E) /m 2 s- 1 = 1 .98 x 10" 1 - 1 .45 x 10~ 3 ■ E + 1 .01 x 10~ 5 ■ E 2 - 1 .83 x 10~ 8 ■ E 3 

D L (E) /mV 1 = 1.33 x 10" 1 +8.87 x 10~ 5 -£+4.40 x 10~ 6 ■ E 2 + 1.10 x 10~ 9 ■ £ 3 

and for the flux coefficients: 

p(E) /m 2 V -1 s _1 = exp [-3.52 - 3.86 x 10" 2 - In £ + 3.26 x 10 _1 /£-(2.57/£) 2 ] 

a*(E) /m _1 = exp [9.67 + 6.46 x 10" 1 -ln£- 4.18 x 10/£] 

£>*(£) /mV 1 = 1 .96 x 10" 1 - 1 .38 x 10" 3 ■ E + 9. 19 x 10~ 6 ■ E 2 - 1 .80 x 10~ 8 ■ E 3 

D* L (E) /mV 1 = 1 .37 x 10" 1 - 4.20 x 10~ 4 ■ E + 6.93 x 10~ 6 ■ E 2 - 1 .56 x 10~ 8 ■ E 3 

ki(E) /m = 6.98 x 10~ 7 + 1.45 x 10~ 2 /(£ + 6.95 x 10) 2 . (A.4) 

Note that the bolsig+ results for argon is only presented from 5 kV/cm to 180 kV/cm. Good 
agreements have been achieved between our flux coefficients and the bolsig+ coefficients in 
electron mobilities, ionization rate. A remarkable difference remains for the transversal diffusion. 
There is also a good agreement between the two sets for the mean electron energies at low field, 
but the difference become to noticeable above 50 kV/cm. 
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